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LONG-TERM  GOALS 

Our  long  term  goal  in  this  and  related  projects  is  to  maximize  the  amount  of  information  regarding 
river  flow  and  morphology  that  can  be  deduced  from  remotely  sensed  information.  In  this  specific 
project,  our  goal  is  to  test  and  improve  the  river  bathymetry  prediction  methods  developed  in  previous 
ONR-funded  work.  This  goal  involves  further  development  on  computational  methods  as  well  as 
exploration  and  refinement  of  methods  for  remotely  measuring  water-surface  elevations  and  water- 
surface  velocities.  In  previous  work,  we  have  been  able  to  demonstrate  that  inversion  of  the 
momentum  equations  and  morphodynamics  modeling  can  be  used  to  detect  errors  and  fill  gaps  in 
various  remotely  sensed  river  bathymetry  data  sets  (i.e.,  such  as  those  measured  using  bathymetric 
LiDAR  and  various  optical  correlation  techniques  using  multi-  and  hyperspectral  scanning,  as  reported 
in  Wright  and  Brock  (2002),  Kinzel  et  al.,  (2007),  Legleiter  and  Roberts  (2005),  and  Legleiter  et  ah, 
(2004)).  In  the  current  project,  we  are  looking  at  this  problem  in  a  more  expansive  sense  by  testing 
measurement  techniques  and  computational  approaches  that  can  be  used  in  the  absence  of  any 
bathymetric  estimate,  as  would  typically  be  the  case  in  rivers  with  high  suspended  sediment  loads 
and/or  deep  depths,  which  limit  optical  techniques.  Eventually,  we  envision  that  these  two  different 
techniques  can  be  rejoined  in  a  manner  that  is  accurate,  efficient  and  capable  of  providing  meaningful 
error  estimates  for  predicted  morphology. 

OBJECTIVES 

The  specific  objectives  for  the  first  year  of  this  two-year  effort  were  to  (1)  modify  our  existing 
inversion  code  to  include  higher-order  effects  neglected  in  our  original  development  (Nelson,  2012), 
and  (2)  to  develop  and  test  methods  for  measuring  water-surface  elevation  and  water-surface  velocities 
in  the  laboratory  and,  if  possible,  at  field  sites.  As  noted  in  our  original  proposal,  we  expected  to 
investigate  two  methods  for  water-surface  elevation  measurements,  one  for  smaller  scale  systems 
(laser-scanning)  and  one  for  larger  systems  (the  AirSWOT  methodology,  see  Rodriguez,  2010).  For 
water-surface  velocity  measurements  we  planned  to  use  the  infrared  videography  methodology 
developed  by  others  (e.g.,  Dugan  and  Piotrowski,  2003;  Dugan  et  al,  2013;  Kinzel  et  al.,  2011).  The 
first  objective  was  relatively  straight-forward,  but  the  second  required  procuring  laser  scanning  and 
infrared  videography  equipment  and  developing  our  capabilities  in  this  area.  Thus,  most  of  our  efforts 
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in  this  first  year  were  directed  at  building  and  testing  an  infrared  videography  system  and  testing 
software  for  deriving  water-surface  velocities  from  that  system.  Our  first-year  work  also  included  a 
smaller  effort  aimed  at  testing  laser-scanning  as  a  method  for  measuring  water-surface  elevation,  and  a 
relatively  minor  programming  element  aimed  at  adding  higher-order  elements  to  the  inversion 
technique  developed  under  our  last  ONR  research  effort. 

APPROACH 


In  our  last  round  of  ONR  research  effort,  we  developed  a  simple  method  for  estimating  depth  based  on 
water-surface  elevation  and  water-surface  (or  vertically  averaged)  velocity.  Instead  of  using  a 
computationally-intensive  inversion  that  solves  for  depth  based  on  a  full  partial  differential  equation  or 
uses  a  data  assimilation  technique  involving  many  flow  solutions,  the  equations  expressing 
conservation  of  mass  and  momentum  were  simplified  based  on  a  few  reasonable  assumptions  to  give 
expressions  that  directly  predict  local  flow  depth  as  a  function  of  water-surface  elevation  and  vertically 
averaged  velocity.  For  example,  from  the  streamwise  component  of  the  equation  expressing 
conservation  of  momentum  for  river  flow,  if  we  assume  the  flow  is  incompressible  and  quasi-steady, 
has  no  lateral  stresses,  and  that  a  simple  drag  coefficient  closure  is  suitable,  we  can  directly  obtain  an 
equation  relating  local  depths  to  local  hydraulic  variables,  as  follows: 


-QV<f/>2  +<v>2  <U> 


g  dE  <u>d<u>  d <u>  <uxv> 

1  -Nds  1  -N  ds  dn  (l -N)R 


where  N  is  a  simple  metric  coefficient,  R  is  the  centerline  radius  of  curvature  of  a  channel-fitted 
centerline,  g  is  gravity,  Cd  is  a  drag  coefficient,  <u>  and  <v>  are  vertically  averaged  velocity 
components  in  the  streamwise  and  cross-stream  directions,  and  s  and  n  are  coordinates  locally  oriented 
in  the  streamwise  and  cross-stream  directions,  respectively. 

To  test  the  expressions  and  investigate  potential  errors,  we  initially  applied  them  to  data  sets  developed 
from  existing  river  models.  Measured  bathymetry  was  used  to  compute  flows  of  a  certain  discharge  in 
river  channels  using  open-source,  public  domain  model  FaSTMECH  (Nelson  and  McDonald,  1996; 
Nelson,  Bennett  and  Wiele,  2003;  McDonald  et  al,  2006)  available  in  the  International  River  Interface 
Cooperative  (iRIC;  www.i-ric.org )  river  modeling  package  developed  by  our  group  and  others.  The 
predicted  water-surface  elevations  and  vertically  averaged  velocities  were  extracted  from  the  model 
results  and  the  depth  was  predicted  from  that  information  using  the  formulation  above.  While  this  is 
clearly  circular,  it  allows  us  to  start  with  data  sets  that  are  constrained  to  satisfy  conservation  of  mass 
and  momentum  and,  as  discussed  in  one  of  our  publications  (Nelson  et  al,  2012),  it  provides  some 
good  physical  guidance  on  issues  with  the  approach  and  the  effect  of  noise  in  the  driving  data. 

We  developed  results  of  this  process  for  two  very  different  river  reaches:  the  Trinity  River  in 
California,  which  is  a  relatively  steep,  shallow,  coarse-bedded  channel,  and  the  Kootenai  River  in 
Idaho,  which  is  a  relatively  flat,  deep,  sand-bedded  reach.  For  each  river,  the  inversion  technique  was 
applied  directly  to  compare  the  computed  and  measured  depths.  The  impacts  of  normally  distributed 
noise  in  the  input  data  sets  (water-surface  elevation  and  velocity)  was  also  explored.  Finally,  the 
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impact  of  using  surface  velocity  instead  of  vertically  averaged  velocity  in  the  inversion  method  was 
investigated.  These  results  are  presented  in  Nelson  et  al.  (2012). 

As  part  of  that  earlier  investigation,  we  were  able  to  show  that  some  terms  neglected  in  the 
development  of  equation  (1)  were  actually  important  in  certain  areas  (notably  lateral  stresses  near 
channel  banks  and  other  regions  of  strong  lateral  shear).  Thus,  as  part  of  our  first  year  effort  on  the 
current  project,  we  developed  a  method  that  brings  these  neglected  terms  in  iteratively.  Thus,  we  now 
solve  Equation  1,  then  use  the  estimated  depths  (and  known  water’s  edge  data)  to  bring  in  the 
neglected  terms  in  a  simple  manner. 

The  next  step  in  this  development  was  to  get  good  enough  laboratory  and  field  data  to  test  the  approach. 
Most  of  our  effort  in  this  first  year  was  spent  on  this.  We  believed  from  previous  work  that  we  could 
get  water-surface  elevations  using  laser  scanners  for  small  scale  field  sites  and  interferometric  radar 
techqniques  (AirSWOT)  for  larger  ones,  but  we  were  much  less  certain  of  our  ability  to  get  water- 
surface  velocities.  We  investigated  a  few  possibilities  and  decided  that  surface  velocity  measured  using 
infrared  videography  (Dugan  and  Piotrowski,  2003;  Dugan  et  al,  2013))  looked  promising  for  remote 
collection  of  the  velocity  field.  The  basic  idea  is  as  follows.  Flowing  water  has  small  temperature 
variations  across  the  water  surface  that  are  coherent  over  timescales  of  seconds  to  minutes.  Two 
examples  are  shown  in  Figure  1 .  These  structures  advect  downstream  at  approximately  the  surface 
velocity  of  the  flow  and  can  be  imaged  with  an  infrared  (IR)  camera.  Our  goal  is  to  create  an  image 
sequence  (a  movie)  of  these  temperature  variations  and  use  a  motion  tracking  algorithm  to  estimate  the 
water-surface  velocity  field  based  on  the  change  in  the  position  of  thermal  structures  between  frames 
of  the  image  sequence. 


Figure  1.  Infrared  images  of  water  surfaces,  a)  an  irrigation  canal  and  b)  a  laboratory  flume. 
Brighter  areas  are  warmer  than  darker  areas  by  <1°  C. 


The  basic  concept  behind  a  variety  of  image  tracking  algorithms  is  that  the  difference  in  the  position  of 
objects  in  two  successive  two  frames  of  a  movie  is  a  displacement  field.  This  concept  is  represented 
mathematically  by  the  Displaced  Frame  Difference  (DFD)  equation: 

I(x,y,t  +  At)~  I(x,y,t)  =  v(x,y)A/  2 
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Where  /  is  the  matrix  of  image  intensity  at  the  pixel  locations  (x,jy)  at  time  t ,  At  is  the  time  elapsed 
between  successive  frames,  and  v(x,y)  =  dl(x,y,t)/dt  is  the  vector  velocity  field  representing  the 
motion  of  objects  in  the  image.  We  explored  two  approaches  to  finding  the  2-dimensional 
velocity  v(x,_y),  particle  image  velocimetry  (PIV)  and  optical  flow  algorithms. 

PIV  is  a  technique  for  imaging  fluid  flows  seeded  with  small  particles  that  move  at  approximately  the 
same  velocity  as  the  fluid  flow.  The  flow  is  illuminated  (typically  with  a  laser)  and  a  series  of  images 
is  recorded.  The  velocity  field  is  estimated  by  finding  the  displacement  that  maximizes  the  correlation 
between  two  successive  frames  of  the  images  sequence.  Typically,  this  accomplished  by  dividing  the 
images  into  m  x  n  pixel  sub-regions  and  finding  a  spatially  averaged  velocity  for  each  sub-region.  See 
Raff  el  et  al.  [2007]  for  a  detailed  discussion  of  PIV  techniques.  We  found  that  PIV  does  not  extract 
uniformly  accurate  velocity  fields  from  the  IR  image  sequences  because  of  the  variations  in  flow 
structure  size. 

Optical  flow  describes  a  broad  class  of  algorithms  that  originate  in  computer  vision  research  for 
tracking  the  motion  of  objects.  Optical  flow  algorithms  use  minimization  schemes  to  find  the  velocity 
field  that  minimizes  the  difference  between  two  frames  on  a  pixel-by-pixel  basis.  Heitz  et  al.  [2010] 
recently  reviewed  optical  flow  algorithms  designed  specifically  for  estimating  fluid  flows.  Corpetti  et 
al.  [2006]  and  Derian  et  al.  [2012]  describe  the  specific  implementations  that  we  have  applied  to  our 
image  sequences.  The  DenseMotion  implementation  described  by  Corpetti  et  al.  [2006]  is  available  as 
a  compiled  Linux  executable  from  http://fluid.irisa.fr/publis/Deliverable2.4.tar.gz.  The  implementation 
of  the  wavelet  technique  described  by  Derian  et  al.  [2012]  is  not  publically  available,  but  Pierre  Derian 
(currently  at  California  State,  Chico)  has  been  collaborating  with  us  and  running  some  analysis.  We 
expect  that  his  software  will  become  available  at  some  point  in  the  future  and  he  has  expressed  interest 
in  further  collaboration. 

The  camera  we  used  is  a  FL1R  SC  8300  HD.  The  camera  has  an  indium  antimonide  (InSb)  sensor  that 
detects  mid-range  IR  light  in  the  3-5  micron  (pm)  band.  The  sensor  size  is  1280  x  720  pixels. 
Temperature  differences  of  25  milli-Kelvins  (0.025  K)  result  in  a  signal-to-noise  ratio  of  1  on  the 
camera  sensor.  The  minimum  resolvable  temperature  difference  is  few  times  this  number  (on  the  order 
of  0.1  °C). 

WORK  COMPLETED 

Based  on  research  efforts  over  the  first  9  months  of  the  first  year  of  our  ONR  grant,  the  following  tasks 
are  complete: 

(1)  We  developed  a  computational  method  for  including  higher-order  terms  in  the  simple  inversion 
method  for  estimating  bathymetry. 

(2)  We  tested  and  confirmed  that  laser-scanners  could  potentially  provide  water-surface  elevation 
data  at  the  required  density  and  accuracy  for  depth  computation  based  on  our  inversion  method. 

(3)  We  purchased  an  infrared  camera  using  ONR  and  USGS  funding  and  developed  a  system  for 
using  it  in  the  laboratory  and  field  to  collect  infrared  videography  for  measuring  water-surface 
velocity. 
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(4)  We  carried  out  several  laboratory  experiments  using  the  infrared  camera  to  measure  water-surface 
velocities 

(5)  We  tested  several  PIV  methods  and  two  optical  motion  methods  for  computing  surface  velocities 
from  infared  videography. 

(6)  We  carried  out  three  field  experiments  to  measure  water-surface  velocity,  two  on  local  canals 
(Croke  Canal  and  Church  Ranch  Canal)  and  a  third  on  the  Platte  River  near  Confluence  Park  in 
Denver. 

RESULTS 

In  this  section,  we  concentrate  on  the  results  of  the  water-surface  velocity  measurement  work.  One  of 
the  first  things  we  learned  from  our  measurement  program  was  that,  although  PIV  could  provide 
velocity  estimates  from  the  infrared  imagery,  it  was  difficult  to  get  uniform  measurement  coverage. 
This  was  because  of  the  importance  of  structure  size  to  the  choice  of  an  optimum  PIV  interrogation 
window.  Thus,  a  PIV  interrogation  window  that  works  for  certain  structure  sizes  fails  for  much  smaller 
ones,  so  for  different  areas  of  the  same  image,  there  is  no  single  window  that  works  well.  This  means 
that  it  is  difficult  to  get  a  complete  coverage  field  of  accurate  velocity  estimates  unless  all  the  observed 
structures  are  similar  in  size  and  the  interrogation  window  can  be  appropriately  tuned.  We  these 
lessons  learned,  we  turned  to  more  general  optical  flow  motion  methods,  as  discussed  here. 

We  recorded  a  series  of  image  sequences  of  the  flow  in  the  small  flume  at  the  U.S.  Geological  Survey 
Geomorphology  and  Sediment  Transport  Laboratory  (GSTL).  The  IR  camera  was  mounted  on  a  rack 
~lm  above  the  surface  of  the  flow  and  oriented  so  that  the  image  plane  was  approximately  parallel  to 
the  surface  of  the  water.  We  assumed  that  the  image  plane  was  exactly  parallel  to  the  water  surface  and 
that  there  is  no  distortion  across  the  view  so  the  scale  (in  pixels  per  cm)  is  uniform  across  the  image. 

In  the  laboratory,  where  water  temperatures  tend  to  be  very  uniform  and  viscous  heating  is  very  low, 
we  obtained  the  best  results  using  a  combination  of  ice  and  aquarium  heaters  to  “seed”  the  flow  with 
temperature  gradients.  We  floated  a  10-pound  block  of  ice  (in  a  bag)  in  the  flume  upstream  of  the 
camera.  Immediately  downstream  of  the  ice  (and  upstream  of  the  camera),  we  attached  two  200  watt 
and  two  300  watt  aquarium  heaters  to  the  side  of  the  flume  and  placed  a  ~15  cm  tall  Plexiglas  block  on 
the  bottom  of  the  flume  to  increase  vertical  motion  in  the  fluid.  One  image  from  the  sequences  is 
reproduced  in  Figure  lb.  A  ruler  placed  on  the  side  of  the  flume  where  it  was  visible  in  the  (un¬ 
cropped)  images  provided  scale  control. 

The  image  sequence  was  recorded  at  a  resolution  of  384  x  510  pixels  at  75  frames  per  second  (fps).  A 
minimal  number  of  frames  (~2)  were  dropped.  The  image  sequence  was  cropped  to  320x320  pixels  (to 
reduce  processing  time)  and  filtered  both  spatially  and  temporally  to  reduce  noise.  We  also  enhanced 
the  contrast  with  an  image  processing  application. 

In  general,  the  wavelet  algorithm  described  by  Derian  et  al.  [2012]  gives  a  better  estimate  of  the 
velocity  field  than  that  of  Corpetti  et  al.  [2006].  This  can  be  seen  in  Figures  2  and  3,  which  show 
fields  of  virtual  “particles”  moving  according  to  the  estimated  velocity  fields.  In  the  wavelet  estimate 
(Figure  2),  “particles”  follow  the  advection  of  the  thermal  structure  more  faithfully  than  in  the 
DenseMotion  estimate  (Figure  3).  This  is  easier  to  see  on  the  videos  on  our  ONR  QUAD  than  it  is 
using  a  still  image,  as  shown  here.  The  velocity  distribution  estimated  by  the  wavelet  algorithm  is 
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shown  in  Figure  4.  The  DenseMotion  estimate  is  shown  in  Figure  5.  No  independent  detailed 
measurement  of  the  flow  velocity  exists,  but  the  visual  correspondence  between  the  motion  of  thermal 
structures  and  the  motion  of  the  “particles”  in  Figure  2  suggests  that  the  wavelet  estimate  is  accurate. 
Both  the  wavelet  and  DenseMotion  algorithms  have  difficulty  estimating  velocity  in  areas  of  low 
image  contrast,  but  the  problem  is  much  worse  in  the  DenseMotion  estimate,  which  results  in 
clustering  of  the  virtual  particles  (Figure  3)  at  high  contrast  boundaries  and  an  overall  lower  estimate 
of  the  mean  velocity. 


Figure  2.  Visualization  of  velocity  field  estimated  with  the  wavelet  algorithm  using  particles 
advected  by  the  predicted  velocity  field.  See  QUAD  for  video  animation. 
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Figure  3.  Visualization  of  velocity  field  estimated  by  DenseMotion  using  particles  advected  by 
velocity  field  estimate.  See  QUAD  for  video  animation. 


Figure  4.  Velocity  distribution  estimated  by  the  wavelet  algorithm.  Cumulative  distribution  function 
(CDF)  is  on  the  left,  probability  density  function  (PDF)  is  on  the  right.  The  red  lines  are  the  best 

fitting  Gaussian  distribution. 
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Figure  5.  Velocity  distribution  estimated  by  DenseMotion.  Cumulative  distribution  function  (CDF) 
is  on  the  left,  probability  density  function  (PDF)  is  on  the  right.  The  red  lines  are  the  best  fitting 

Gaussian  distribution. 


In  addition  to  many  similar  laboratory  sequences,  we  also  deployed  the  camera  system  in  three  field 
settings.  Figure  6  shows  an  image  from  a  sequence  taken  from  a  bridge  over  the  South  Platte  River  at 
Confluence  Park  in  Denver  in  the  early  morning  of  July  31,2013.  The  visualization  of  the  wavelet 
velocity  field  estimate  using  advected  particles  is  shown  in  Figure  7  (a  video  of  this  interesting  flow 
field  is  also  shown  in  our  ONR  QUAD  sheet). 


8 


Figure  6.  South  Platte  River  recorded  from  a  pedestrian  bridge  in  Confluence  Park,  Denver  with  a 

visualization  of  the  wavelet  estimate  of  the  flow  field. 


We  encountered  a  variety  of  issues  during  this  work.  First,  infrared  radiation  reflected  off  of  the  water 
surface  is  a  significant  problem.  On  sunny  days,  the  reflection  of  the  sky  off  the  water  surface  makes 
the  thermal  structure  much  more  difficult  to  image.  This  is  less  of  a  problem  on  overcast  days.  The 
irrigation  canal  image  shown  in  Figure  1  was  shot  at  mid-day  on  an  overcast  day  in  February.  The 
image  sequence  shown  in  Figure  7  and  discussed  in  the  next  section  was  recorded  in  the  mid-afternoon 
on  an  overcast  summer  day  (8/7/13). 
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IR  glare  can  be  a  problem  even  at  night.  Glare  is  most  evident  if  there  is  vegetation  along  the  banks 
and  the  channel  is  narrow.  We  believe  this  is  due  to  infrared  radiation  from  the  vegetation.  Because 
light  reflected  off  of  non-metallic  reflective  surfaces  becomes  polarized,  we  believe  some  of  the  issue 
with  reflected  radiation  can  be  mitigated,  but  we  have  not  yet  tried.  Visual  light  photographers  use 
polarizing  fdters  to  reduce  glare.  It  seems  reasonable  to  expect  that  an  IR  polarizing  filter  might  be 
similarly  useful.  Infrared  polarizers  are  available  for  a  few  thousand  dollars. 


Figure  7.  Irrigation  canal  recorded  on  the  overcast  afternoon  of  8/7/13. 

In  addition  to  exploring  these  techniques  for  measuring  detailed  surface  velocities,  some  preliminary 
work  on  measuring  field  water-surface  elevations  with  laser  scanners  was  completed.  This  was  done 
because,  even  though  we  had  accurate  methods  for  measuring  detailed  water-surface  elevation  in  the 
lab,  the  techniques  are  not  amenable  to  application  over  large  field  areas;  using  a  laser  scanner  is 
potentially  a  good  field  method,  although  we  were  unable  to  find  previous  work  where  this  method  had 
been  applied  and  tested  with  other  measurements.  We  carried  out  several  field  data  acquisitions  and 
were  able  to  show  that  the  scanner  can  measure  surface  elevations  very  quickly  over  an  entire  reach,  as 
shown  in  Figure  8. 
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Figure  8.  Laser  scanner  image  of  water-surface  topography  at  the  confluence  of  the  Platte  River 
and  an  irrigation  canal.  Comparisons  to  ground-based  surveying  suggests  we  can  resolve  the  water 

surface  to  within  l -2mm,  less  with  averaging. 

Finally,  as  part  of  the  same  effort  to  develop  field  techniques  for  measuring  water-surface  elevation, 
we  spoke  with  researchers  on  the  NASA  AirSWOT  program  and  were  instrumental  in  planning  a 
combined  effort  for  data  collection  on  the  Sacramento  River  this  Fall.  AirSWOT  surface-elevation  data 
will  be  collected  and  USGS  will  work  on  collecting  both  truth  data  for  comparison  and  bathymetry  and 
velocity.  We  hope  to  include  some  component  of  infrared  videography  for  water-surface  velocities,  but 
this  is  uncertain  at  this  point. 

IMPACT/APPLICATIONS 

The  research  work  carried  out  under  this  project  is  also  playing  a  significant  role  in  USGS 
developments  for  practical  applications  of  non-contact  gauging.  In  addition,  we  have  added  the 
inversion  model  code  to  the  iRIC  modeling  interface,  so  it  will  be  freely  available  to  the  public  in  the 
next  major  release,  expected  May  2014. 

RELATED  PROJECTS 

We  continue  to  tie  this  work  to  other  research  projects,  including  work  on  the  Kootenai,  Platte,  and 
Russian  Rivers.  In  addition,  we  are  working  closely  with  NASA  and  personnel  at  the  USGS  California 
Water  Science  Center  on  the  application  of  the  AirSWOT  water-surface  topography  testing  program 
on  the  Sacramento  and  hope  to  have  access  to  those  data  sets  sometime  later  this  year  or  early  next. 
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